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ABSTRACT 


Using  an  equation  of  state  of  air  in  completely  tabular  form, 
a  one  dimensional,  spherically  symmetric  blast  wave  calculation  has 
been  numerically  carried  out.  An  IBM  704  computer  was  utilized  for 
the  calculation.  The  Rankine-Hugoniot  conditions  at  the  shock  front  and 
the  isentropic  changes  of  the  shocked  fluid  were  determined  by  iterative 
methods.  The  numerical  methods  employed  are  discussed  in  some  detail, 

V 

as  are  the  details  of  the  equation  of  state.  The  initial  starting  con- 
*■  ditions  for  this  numerical  integration  were  those  of  Problem  M.  The 

numerical  results  are  presented  in  graphical  form.  Comparison  to  Problem 
M  is  also  displayed  graphically.  The  comparison  to  Problem  M  is  good, 
the  small  differences  being  attributed  to  that  of  equation  of  state  and 
finite  differencing  methods. 
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1.  INTRODUCTION 


This  work  was  undertaken  in  order  to  develop  a  sharp  shock 
blast  wave  code.  As  a  means  of  checking  the  present  coding,  Problem  M,^ 
a  blast  wave  calculation  for  an  energy  release  of  13«5  KT  carried  out 
at  Los  Alamos  around  19^5,  was  chosen  as  the  comparison  problem.  The 
initial  conditions  of  Problem  M  were  used  as  the  starting  data  for  this 
code,  thereby  permitting  each  calculation  cycle  to  be  compared  to  those 
of  Problem  M.  Such  a  comparison  serves  as  a  guide  to  determine  if  the 
present  code  is  functioning  properly,  and  if  it  is,  serves  also  as  a 
check  on  Problem  M,  for  since  the  results  of  that  problem  have  been  used 
extensively,  an  independent  comparative  blast  calculation  is  desirable. 
The  results  from  this  problem  differ  somewhat  from  those  of  Problem  M, 
since  the  calculations  for  the  latter  were  carried  out  on  slow,  semi¬ 
automatic  equipment,  with  a  substantial  portion  of  M  done  by  hand  on 
desk  calculators.  Thus  approximations  consistent  with  realistic  calcu¬ 
lation  times  were  necessary.  In  addition,  a  new  equation  of  state  for 
2  5  k 

air  ’  *  was  used.  Improvements  in  the  numerical  methods  of  Problem  M 
were  possible  using  the  present  day  high  speed,  completely  automatic 
computing  machines. 

The  new  equation  of  state  data  were  in  tabular  form.  Several 
exploratory  attempts  to  find  an  analytic  fit  to  these  data  indicated  that 
a  complex  system  of  fits  would  undoubtedly  be  necessary.  In  general, 
since  most  analytic  fits  to  equation  of  state  data  are  poor  except  in 


limited  regions,  and  since  the  range  of  thermodynamic  variables  in  blast 
problems  is  quite  extreme,  it  was  thought  that  if  this  new  equation  of 
state  were  used  completely  in  tabular  form,  it  would  be  more  useful 
and  accurate  than  an  extremely  complex  system  of  fits.  Given  the  rela¬ 
tive  specific  volume  and  the  specific  internal  energy,  this  tabular 
equation  of  state  determines  the  pressure  by  a  double  interpolation  of 
the  tabular  entries,  or  symbolically,  P  =  P(v/Vq,E).  The  details  of 
the  equation  of  state  are  discussed  in  Section  9* 

The  method  used  to  determine  the  isen tropic  changes  within 
the  fluid  of  the  shocked  sphere  is  given  in  Sections  2  and  3. 

The  Rankine-Hugoniot  conditions  at  the  shock  boundary  are 
determined  by  an  iterative  method  which  is  discussed  in  Section  4. 

The  addition  of  new  fluid  elements  to  the  calculation  is  re¬ 
quired  as  the  shock  discontinuity  progresses  further  into  the  undis¬ 
turbed  fluid.  This  addition  is  considerably  complicated  by  utilization 
of  a  "sharp  shock"  calculation'  >z>’  rather  than  a  "smeared  shock"  calcu¬ 
lation.^®’^  The  addition  procedure  is  explained  in  Section  5,  and 
other  fluid  element  adding  schemes  which  were  tried  are  described  in 
Section  6. 

The  initial  conditions  chosen  for  this  numerical  integration 
represent  the  hydrodynamic  state  of  an  energy  release  of  13*2  KT  at  a 
time  12  milliseconds  after  detonation.  There  is  a  difference  in  total 
energy  between  this  problem  and  that  of  Problem  M  because  of  the  new 
equation  of  state  used.  Comparison  values  from  the  present  problem  and 
from  Problem  M  are  displayed  graphically  in  Section  11. 


2.  THE  BASIC  EQUATIONS  AND  BOUNDARY  CONDITIONS 


We  consider  the  region  of  a  fluid  bounded  by  a  shock  front. 
Within  this  bounded  region  the  hydrodynamical  state  of  the  fluid  is 
described  by  four  equations:  the  equation  of *  continuity,  the  equation 
of  motion,  the  conservation  of  energy  in  the  form  of  the  first  law  of 
thermodynamics,  and  the  equation  of  state.  Spherical  symmetry  is  assumed 
and  the  Lagrangian  form  of  the  equation  of  motion  is  used.  In  the  Lagran 
gian  formulation  we  are  concerned  with  the  properties  of  fluid  elements, 
which  are  tagged  or  identified  by  their  initial  positions,  and  we  follow 
these  elements  along  in  space  and  time  observing  the  changes  of  their 
fluid  properties.  Thus,  the  radius,  pressure,  density,  velocity,  in¬ 
ternal  energy,  and  acceleration  of  each  fluid  element  are  functions  of 
the  element's  initial  position  and  the  time,  or  R(r,t),  P(r,t),  p(r,t), 
v(r,t),  E(r,t),  and  a(r,t).  R  is  referred  to  as  the  Eulerian  radius 
and  r  as  the  Lagrangian  radius,  that  is,  R  is  the  physical  position  of 
a  fluid  element  at  time  t,  and  r  is  that  element's  initial  position. 

The  equation  of  continuity  or  conservation  of  mass  states 
that  the  initial  mass  of  a  specific  fluid  element  remains  in  that 
element  and  there  is  no  transport  of  mass  from  one  element  to  another, 
or 


fo 

P 


(2.1) 


where  V  is  the  specific  volume  or  the  reciprocal  of  the  density. 


From  Newton’s  second  law  the  volume  force  is  given  by 


F 


V 


-grad  P 


and  for  spherical  symmetry. 


dv  _  „i  ap 
5t  p 

Combining  this  with  (2.l),  we  have 

dv  =  1  ]£  3P  =  3R 2 

^  =  "  P0  r2  *  "  "  p0  Sr5 


(2.2) 


(2.2’) 


We  require  the  changes  of  the  fluid  within  the  bounded  volume  to  be 
isentropic,  or 


dE  =  -PdV. 


(2.3) 


The  equation  of  state  is  in  tabular  form.  The  table  is  arranged  so 
that  for  a  given  relative  specific  volume  and  specific  internal  energy 
the  pressure  is  determined  by  a  double  interpolation  scheme  or,  symboli¬ 
cally, 

P=p[v/v0,e].  (2-^) 

See  Section  9  for  a  more  detailed  discussion  of  the  equation  of  state. 

The  above  equations  are  valid  behind  the  shock.  At  the  shock, 
however,  the  equations  of  Rankine  and  Hugoniot  are  to  be  used.  They  are 
as  follows: 
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p(u  -  l)  =  P0I  ,  (2.5) 

*  -  P0  -  P0£u  ,  (2.6) 

Eg  -  Eq  =  £(*  +  PQ)  (l/p0  -  l/pg).  (2.7) 

The  material  velocity  of  the  unshocked  region,  u^,  is  assumed  to  be  zero. 
i  is  the  shock  velocity,  u,  the  material  velocity  of  the  medium  behind 
the  shock,  and  \|r,  the  shock  pressure. 

There  are  two  boundary  conditions,  one  at  the  center  of  the  dis¬ 
turbance,  the  other  at  the  shock.  They  are  as  follows: 

R(0,t)  =  0  (for  all  t)  (2.8) 

and 

R(!,t)  =  £  (for  all  t),  (2.9) 

where  £  is  the  Lagrangian  radius  of  the  shock  front  and  is  a  function 
of  time  or  £  =  £(t).  Equations  (2.5),  (2.6),  (2.7),  and  (2.9)  are  not 
independent , as  can  be  readily  verified. 

3.  INTEGRATION  OF  THE  BASIC  EQUATIONS 

The  integration  of  the  basic  equations  is  carried  out  numeri¬ 
cally  in  a  stepwise  manner.  The  differential  equations  are  approximated 
by  finite  difference  equations.  In  order  to  set  up  these  difference 
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equations  the  hounded  region  is  divided  into  finite  regions  or  fluid  ele¬ 
ments  of  concentric  spherical  shells  each  of  equal  Lagrangian  radial  width. 


The  radii  bounding  the  elements  are  numbered  consecutively  outwards  from 
the  center,  0,  1,  2,  3,  4,  •••,  i-1,  i,  i+1, • • • •  pressure,  density, 

and  internal  energy  of  the  entire  element  are  considered  to  be  that  of  the 
centroid  of  each  fluid  element  but  the  acceleration  and  velocity  are  asso¬ 
ciated  with  the  bounding  radii.  Let  us  consider  the  equations  for  a  typical 
fluid  element  bounded  by  radii  ^  and  r^  at  time  t  =  n-1  and  integrate 
these  equations  over  a  time  interval  At.  We  let  the  subscripts  indicate 
the  Lagrangian  radii  where  the  quantity  is  assumed  located  and  the  super¬ 


script, 


the  time. 

n-1 

a.. 


Equation  (2.2')  then  becomes 

fife'1)2  ^3  -  pr-*  v-4  - 

■  "o  3  7T~  ~  At 

p0  Ar.,1  +  Ar.  i 

1+2  1-2 


(3-1) 


where  Ar?  x  =  (r?  ,  -  r?),  and  where  an  average  of  the  two  bounding 
i+2  i+l  i 

Ar^'s  is  utilized  for  computational  convenience.  Solving  for  the  velocity. 


we  have 


vn-i  .  yn-3/2  +  an-l 
i  i  i 


g  -  *rx 

At 


and  for  the  radius. 


Ri  =  Ri_1  +  yn"2  At  * 


(3-2) 


(3-3) 


Now,  having  the  new  outer  radius  of  the  element  at  t  =  n  and  assuming 
that  we  have  previously,  in  a  similar  manner,  determined  the  inner  radius 
at  t  =  n,  we  are  in  a  position  to  determine  the  volume  of  the  element  at 
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* 


a 


time  t  =  n. 


HU  ■  ^7  [(Rif  -  [Rhf] 


1"2 


3 

Ar.  i 
!~2 


(Ri-Ri-i)[(Ri)2  +  Ri-i(Ri  tRli)J 


(3.4) 


Equations  (2.3)  and  (2.4)  in  difference  form  become 


E?  i 

1-2 


=  ^4-^2^-  L(V/V0)^ 


(3-5) 


and 

P  =  p(v/v0,  E )  .  (3.6) 

Assuming  that  we  know  all  the  above  quantities  except  pt?  _i  and  E1}  i 

1-2  1-2 

and  have  determined  the  volume  of  this  element  at  time  t  =  n-1,  we  can 

solve  Equations  (3*5)  and  (3*6)  simultaneously  by  iteration. 

We  have  demonstrated  how  a  typical  element  is  integrated 

from  time  t  =  n-1  to  t  =  n;  this  process  is  carried  out  for  each  element 

starting  from  the  center  and  working  toward  the  shock  front.  The  last 

3 

element,  however,  needs  to  be  treated  differently,  as  there  is  no  Ar 
conveniently  available.  The  acceleration  is  formed  as  follows: 


(3-7) 
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The  subscript  I  indicates  the  element  adjacent  to  the  shock.  After 
carrying  out  the  above  scheme  we  have  all  the  information  required  within 
the  bounded  region.  We  now  determine  the  shock  quantities. 


k.  DETERMINATION  OF  THE  SHOCK  QUANTITIES ,  OR  THE  SHOCK  FITTING 

The  method  used  to  determine  the  shock  conditions  is  discussed 
in  IA-1148  (Ref.  5)  and  is  an  iterative  scheme.  We  briefly  outline  it 
here  as  follows:  One  guesses  a  shock  volume  =  V^/VQ  (the  recipro¬ 

cal  of  the  shock  compression) .  From  the  simultaneous  solution  of 
Equations  (2.U)  and  (2.7)  (an  iterative  process),  the  values  of  \|r,  E^, 
and  =  -S  of  the  Hugoniot  are  determined.  The  shock  velocity  is 

dVi 

=A/Vq  w",  (k.i) 


where  (w11)^ 
t  =  n  by 


The  position  of  the  shock  is  determined  at  time 


.n 


r1  ♦  f  (in 


+ 


(U.2) 


Using  the  relation  for  the  total  time  derivative  of  the  shock  volume  in 
terms  of  the  above  determined  quantities,  we  have 


and  in  difference  form, 


Ik 


From  the  above  and  the  value  at  n-1  ve  can  determine  an  average  change 

in  the  shock  volume  in  the  time  interval  At.  Adding  this  to  \  we 

have  a  calculated  volume  V?,  or 

c 


AV*^  =  f  + 


(4.4) 


ev“  =  v^'1  +  av^"2  .  (4.5) 

Comparison  of  V?  and  V?  determines  whether  the  iteration  is  com- 
g  £  c  5 

pleted.  If  it  is  not,  a  new  guess  is  carried  out  based  on  previous  guesses. 
When  the  "shock  fitting"  is  completed,  the  integration  of  all  the  equations 
from  time  t  =  n-1  to  time  t  =  n  has  been  completed.  Integration  from 
time  t  =  n  to  time  t  =  n+1,  etc.,  is  carried  out  by  repeating  the  steps 
of  Sections  3  and  4.  However,  when  the  shock  position  |  becomes  greater 
than  r^  +  Ar  =  ri+i>  a  new  fluid  element  must  be  added  to  the  calculation. 


5.  ADDING  A  NEW  FLUID  ELEMENT 


In  order  to  add  a  fluid  element,  or  mass  point,  one  needs  to 
obtain  values  of  the  pressure,  volume,  internal  energy,  velocity,  and 
the  radius  at  the  appropriate  Lagrangian  radii  and  times,  i.e.,  we  need 


Since  the  shock  fitting  is  dependent 
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upon  the  pressure  and  volume  gradients  behind  the  shock  front  [cf*  (M)], 

the  new  fluid  element  values  and  the  shock  quantities  are  determined 

simultaneously  by  an  iteration  scheme.  Thus  during  the  course  of  a  shock 

fitting  when  the  results  of  Equation  (k.2)  indicate  that  |n>  Tj  +  Ar  =  rI+1, 

the  normal  shock  fitting  scheme  is  interrupted  and  the  above  mentioned 

values  of  the  new  fluid  element  are  created  in  a  method  consistent  with 

V..  (This  method  is  discussed  in  more  detail  later  in  this  section.) 
g  I 

The  shock  fitting  is  then  continued  at  {k.j),  using  now,  however,  the 

pressure  and  volume  just  created.  If  the  comparison  of  and  ^  is 

not  satisfactory,  a  new  V,  is  made  and  the  mass  point  adding,  shock 

g  5 

fitting  iteration  is  continued.  The  creation  of  the  new  element  values  is 
carried  out  on  each  of  these  iteration  cycles.  This  iteration  scheme  is 


continued  until  V.  and  V,  agree  to  the  desired  accuracy.  The  inter- 
g  i  c  | 

ruption  of  the  shock  fitting  scheme  is  made  only  on  the  intergration 


cycle  when  +  Ar  =  ri+q‘  After  the  new  fluid  element  is  added, 

the  interruption  of  the  shock  fitting  scheme  stops  and  the  calculation 


continues  in  its  normal  manner,  except  that  there  is  one  more  mass  point, 


until  it  is  again  determined  by  (k.2)  that  another  element  needs  to  be 


added. 


The  quantities  of  interest  for  the  new  fluid  element  are 
created  by  the  following  method.  First,  the  time  at  which  the  shock 


position  equals  rj  +  Ar  =  rJ+1  is  determined. 
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Let 


.n-1 


*1  = 


1+1 


„n 
|  -  r 


In-l 


and 


*2  = 


1+1 


P 


In  general,  since  £  At,  a  weighted  average  of  t^  and  tg  is  used 

n 

to  describe  the  time  of  crossing.  This  time  becomes  t*  =  t  +  T^,  where 


h  -s  (2At  ‘  \  '  V>-  f5-1* 

We  also  designate  Tg  such  that  Tg  =  At  -  T^.  For  convenience  we  indi¬ 
cate  by  I  the  fluid  element  being  added.  The  quantities  of  interest 
for  the  new  mass  point  at  the  time  t  =  t*  axe  now  determined.  The  radius 
is 

R*  =  rT  =  t*  .  (5-2) 

The  shock  volume  at  this  time,  t  =  t*,  is  determined  by 


V*  =  —  (T  Vn-'*'  +  T  V11). 
v£  At  V12  |  X1 V 


(5-3) 


Using  V*,  we  next  solve  Equations  (2.7)  and  ( 2.4 )  for  E*  and  \| f*.  The 
inner  bounding  radius  is  determined  at  time  t  =  t*  by 


R*  =  Rn  -  v11”2  T  , 
1-1  v  1-1  "l*2, 


(5.4) 


and  the  volume  at  time  t  =  t*  is 


V*  1  =  , 

2  *£  ±. 
1-2 


£7  (r!  ■  E!-i)  [4  -  Ef-1  ( R!-i +  ri)j  •  <5-5) 
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The  internal  energy  at  time  t  =  t*  is  obtained  by  a  three  point  interpola¬ 
tion,  or 

E*_i  =  0.533  e*  +  0.667  e*_3/2  -  0.200  e*_5/2,  (5-6) 


where 


and  E|_5/2 


are  determined  by  interpolation  from  values  at 

times  t  =  n  and  t  =  n-1.  P*  1  is  obtained  through  the  equation  of  state 

1-2 

(2.W  since  V*  x  and  E*  1  are  now  known.  The  velocity  of  the  new  mass 
I-'2  I~2 

point  at  t  =  t*  is  the  material  velocity  of  the  shock,  or 


VI 


-  [vo(**  -  p0)  (p  -  vt)r  ■ 


(5.7) 


We  now  know  all  the  necessary  quantities  of  the  new  mass  point  at  t  -  t*; 

however,  in  order  that  they  will  fit  into  the  general  differencing  scheme 

it  is  necessary  to  advance  them  to  their  proper  times.  First  we  determine 

Ra,  va  2,  and  1.. 

11  1"2 

^  =  r"'1  ♦  at,  (5.8) 

where 

Rj”1  =  rT  -  v*^  +  (5-9) 


1  a* 

vn~2  =  v*  +  —  (T 
I  12  v  2 


-  V' 


(5.10) 


and 


r 


al  = 


Pf  1 
1-2 


Ar 
}0  2 


(5*11) 
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The  volume  of  the  new  element  is  then  at  t  =  t 


n 


V?  i 

1-2 


ArJ  i 

J-“2 


{(Rx 


4-l)  [(Klf 


+  RI-l(RI +  RI-l)]) 


(5.12) 


Since  we  have  P*  .1,  E*  .1,  V*  .1,  and  jl>  we  can  solve  Equations  (5*5)  and 

1-2  1*2  -L~2  1”2 

(3.6)  for  the  values  of  i  and  i,  that  is,  we  require  the  newly 

I"  2  1’2 

formed  element  to  change  isentropically  from  t  =  t*  to  t  =  tn.  We  now 
have  all  the  necessary  quantities  for  the  new  point. 


6.  OTHER  ELEMENT  ADDING  SCHEMES 

Several  other  schemes  were  tried  in  order  to  add  the  new  fluid 
element;  however,  that  of  Section  'y  proved  to  he  the  most  satisfactory  in 
regards  to  smoothness  of  the  P,  V,  E,  v,  and  R  profiles  at  the  time  of 
adding,  and  particularly  in  regards  to  the  conservation  of  total  energy. 

The  other  methods  are  now  briefly  discussed.  The  determination 
of  t*  is  the  same  for  all  the  methods  tried  and  is  that  described  by 
Equation  (5*1). 

■Jfc  ^  1 

a.  Use  of  the  Shock  Quantites  at  Times  t  and  t 

n  n-~ 

Equations  (5.8)  -  (5.12)  were  used  to  determine  R  ,  v^  2,  and 
with,  however,  the  exception  that  was  used  for  a*.  Having  the 

*  *-i 

shock  quantities  at  the  present  *  time  t  and  at  the  previous  *  time  t  , 
the  assumption  is  made  that  the  internal  energy  of  the  new  element  at  time 
t  =  tn  is  the  average  of  these  shock  internal  energies.  Then  P^  _i  is 
determined  by  the  equation  of  state  (2.4). 
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This  method  produced  relatively  smooth  profiles  at  the  time  of 
adding;  however,  the  conservation  of  the  total  energy  of  the  problem  was 
not  so  good  as  that  of  the  method  of  Section  5- 


b.  Original  point  on  the  Hugoniot 


n  n — 

Using  the  scheme  of  6a  to  determine  Rj,  Vj  2,  and  Vj_.i,  a 
point  on  the  Hugoniot  between  time  t*  and  time  t**1  is  chosen  as  repre¬ 
sentative  of  the  entire  new  element.  We  follow  along  the  adiabat  or  ex¬ 
pand  the  new  element  isentropically  to  the  volume  v”  j.  In  this  manner 

I~2 

we  determine  the  pressure  P11  ^  and  the  internal  energy  E”  i  at  time 

1-2  "2 

t  s  tn.  The  representative  point  chosen  was  the  point  corresponding  to 
the  shock  quantities  as  the  shock  crossed  the  centroid  of  the  new  element. 
The  results  proved  to  be  poor,  and  changing  the  representative  point 
brought  little  improvement.  This  scheme,  though  physically  real,  is  not 
consistent  with  the  differencing  scheme.  This  is  clear  if  one  considers 


the  region  between  the  shock  and  the  last  mass  point.  This  region  is  not 
carried  (in  this  problem)  in  the  regular  advancing  scheme  until  the  new 
point  is  added,  i.e.,  the  differencing  scheme  of  this  small  region  con¬ 
siders  that  only  the  region  possesses  mass.  The  energy,  pressure,  etc., 
are  ignored  until  the  new  mass  point  is  added,  and  only  then  are  the 
isentropic  changes  carried.  Thus  this  method  is  not  appropriate  for  our 


present  differencing  scheme. 

Carrying  the  region  between  the  shock  and  the  last  mass  point 
as  a  regular  mass  point,  that  is,  considering  isentropic  changes,  would 
be  advantageous;  however,  the  special  calculations  involved  make  it  cum 


bersome  to  handle  such  a  calculation  for  the  numerical  integration. 
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c.  v*  Variation 

Since  it  can  be  argued  that  one  really  needs  an  average  velocity 
for  the  differencing  scheme  and  not  the  bounding  value,  v*  was  varied  in 
schemes  a  and  b  by  utilizing  various  space  weighting  factors.  Although 
this  method  was  effective  in  changing  the  newly  added  values  to  adjust 
smoothly  in  specific  cases,  it  was  not  possible  to  find  a  weighting 
scheme  which  would  be  satisfactory  for  all  cases. 

d.  Gradient  as  Three  Point  Fits 

In  the  shock  fitting  scheme  of  Section  4,  depends  on 


and 


(see  4.3)* 


Equation  (4.3')  uses  a  two  point  difference  as 


the  approximation  to  the  derivative.  A  three  point  fit  was  carried  out 
for  both  V(r)  and  P(r)  and  the  derivative  evaluated  at  r  =  |  determined. 
These  values  were  then  substituted  in  (4.3)  a nd  the  shock  fitting  scheme 
used  this  value  of  V^.  This  scheme  functioned  nicely  on  regular  cycles; 
however,  on  the  add  mass  point  cycles  the  values  of  V,  and  V„  would 
agree  to  only  about  5  per  cent.  The  pressure  gradient  used  in  the  last 
mass  point  acceleration  term  was  also  determined  by  a  three  point  fit  of 
P(r);  i.e.,  from  the  fit  was  used  in  Equation  (3*7)*  This  too  had 

no  appreciable  effect  and  did  not  help  the  poor  convergence  described 
above . 


e.  No  Coupling  in  the  Point  Add  Scheme 

In  this  scheme  the  shock  properties  and  the  new  point  quantities 
were  not  determined  simultaneously.  The  new  point  quantities  were  deter¬ 
mined  after  the  shock  fitting;  however,  the  same  methods  outlined  in 
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Section  5  for  the  point  adding  were  still  used.  This  method  though  good 
was  not  so  good  as  that  of  the  simultaneously  determined  method  of 
Section  5. 

f.  Variation  of  the  Time  Interval  At 

It  was  found  that  even  though  the  stability  conditions  were 
good,  the  new  mass  point  data  could  be  improved  if  the  time  interval  At 
was  reduced. 

g.  Combinations  of  the  Above  Schemes 

Several  combinations  of  the  above  schemes  were  tried.  For 
example,  schemes  a  and  e  were  tried,  also  schemes  a,  c,  and  e;  or  the 
scheme  of  Section  5  with  those  of  e  and/or  c,  etc. 

After  trying  a  number  of  these  combinations,  the  method  of 
Section  5  proved  to  be  the  best. 

7.  STABILITY  CHECKS 

The  Courant  stability  condition  was  employed  in  order  to  en¬ 
sure  a  stable  differencing  scheme.  A  stability  number  was  calculated 
for  each  mass  point  every  integration  cycle.  If  the  stability  number 
exceeded  a  particular  value.  At  of  the  next  calculating  cycles  was 
reduced  by  a  factor  of  2.  Similarly,  if  the  stability  number  remained 
below  another  particular  value  for  several  calculation  cycles,  the  At 
of  the  next  calculation  cycles  could  be  increased  by  a  factor  of  2.  The 
problem  was  automatically  monitored  for  the  stability  time  halving,  but 
that  for  the  time  doubling  was  carried  out  visually  by  the  operator.  The 
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approximate  difference  form  of  the  stability  relation  used  was 


"  e“  -  r" 
1-1 


(7.1) 


In  (7.1),  the  sound  velocity  was  approximated  by  that  of  an  ideal  gas 
with  7  =  1.4.  The  maximum  stability  number  and  its  corresponding  mass 
point  number  were  stored  each  integration  cycle  and  printed  as  so  desired. 


8.  ENERGY  CALCUIATIONS 


The  total  energy  of  the  problem  was  calculated  every  integration 
cycle;  however,  the  energy  calculated  on  cycle  n  was  the  energy  for  cycle 
n-1,  the  reason  being  that  the  velocity  is  centered  in  time  at  n-^r. 

The  kinetic  energy  of  each  mass  point  was  computed  as  follows: 


=  ?M  i 
l-2 


(8.1) 


where  M.  1  =  4/3«p„Ar^  1  (the  mass  of  the  element)  and  where 
1*2  U  1“2 


-1 

2. 

1-2 


with 


n-1 

v. 

1 


n-i  n-3/2 

V  +  V 

i _ i 

2 


The  internal  energy  of  a  mass  point  was  calculated  simply  as 


if1.'}  =  E?"!  M.  1. 
1“2  1~2  1“2 


(8.2) 
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The  energy  within  the  small  region  between  the  shock  and  the  last  mass 
point  is  obtained  by  interpolation.  The  kinetic  energy  is 


(k.e.)"_J  =  m{_i. 


(8.3) 


where 


(-S) 


>  ta 

I  "  2 


n--g-  n-3/2 

_  v„  +  v_ 

n-1  I  I _ 


and  where 


Ci  "  l5""1  '  rx)[(t'1"1)  +  ri(  t“'1  +  ri)]‘ 


The  internal  energy  is 


uj'i  »  Eg'i  M  i, 

|-2  5“2  5-2 


(8.4) 


where 


Et"i  =  e*-1  + 
t~2  I 


1  -  En~i 

_ | _ I  “"2 

n-l  n-1 
-1  RI  +  RI-1 


.  i”-1  .  a"’1 

n-1  I 


The  total  kinetic  energy  is  then  the  sum  of  all  these  elemental 


kinetic  energies,  or 


(K.E. )n"*1  =  (K.E.)-I  +  (K.E.)-i  . 


(8.5) 
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Summing  over  all  the  elemental  internal  energies  also  gives  the  total  in¬ 
ternal  energy;  however,  the  ambient  internal  energy  of  the  medium  must 
first  be  subtracted. 

(u)"-1  =  +  £  <£*  -  V3  VoK6”'1)5  -  4  <8-6) 

The  total  energy  of  the  problem  is  then 

£  =  (K.E.)n_1  +  (U)n_1.  (8.7) 

Figures  13,  14,  and  15  show  the  total  energy,  the  total  internal 
energy,  and  the  total  kinetic  energy,  respectively,  plotted  as  functions 
of  the  time . 


9.  THE  EQUATION  OF  STATE  OF  AIR 

The  equation  of  state  is  in  tabular  form.  The  pressure  is 
determined  as  a  function  of  the  specific  energy  and  the  relative  specific 
volume  through  the  table  searching,  double -quadratic,  interpolating  sub¬ 
routine  1-B  and  a  table  of  thermodynamic  values.  This  table  is  composed 
for  specific  values  of  log10(v/vo)  and  corresponding  to  each  is  a  set  of 
values  of  PV  and  E.  There  are  15  log10(v/VQ)  values  ranging  in  0.2  inter¬ 
vals  from  -1.2  to  +1.6.  For  each  log10(v/v0)  there  are  6l  sets  of  PV, 

E  values.  Each  set  of  PV,  E  corresponds  to  a  specific  temperature,  thus 
6l  temperatures  are  represented.  These  temperatures  range  in  value  from 
200°K  to  2,512,000°K. 
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The  basic  data  sources  for  the  table  were  References  2,  3>  b, 
and  10.  The  temperature  range  is  divided  into  three  parts:  low,  inter¬ 
mediate,  and  high.  Reference  2  was  the  source  for  values  in  the  low 
temperature  range  (200°K  -  1900°K).  As  these  data  were  in  quite  a 
different  form  from  that  desired,  it  was  necessary  to  interpolate  the 
data  and  rearrange  it  to  conform  with  the  above  chosen  tabular  form. 
Reference  3  was  the  source  for  the  intermediate  range  values  (2000°K  - 
15,000°K).  Values  in  the  high  temperature  range  (15>950°K  -  2,512,000°k) 
were  taken  from  Reference  b.  Here  again,  the  data  were  not  in  the  de- 
sired  form  and  interpolation  was  required.  The  contribution  of  radiant 
energy  to  these  thermodynamic  values  was  added  using  the  relation 


E 


r 


JL  5_JlL  3!^  _ 

15  *  (ho)3  p  "  p 


where  a  is  the  radiation  density  constant 


(9.1) 


(PV)r  =  l/3Er.  (9.2) 

After  these  processes  were  carried  out  for  the  high  temperature  range. 
Reference  10  became  available.  Comparison  of  the  values  in  Reference  10 
with  the  above  interpolated  ones  proved  quite  good,  and  the  high  tempera¬ 
ture  range  values  were  therefore  not  changed.  As  each  of  these  references 
assumes  somewhat  different  models  for  the  basic  constituents  of  the  air 
gas,  it  is  not  surprising  to  find  that  in  the  regions  of  overlapping 
there  was  some  disagreement.  Although  the  disagreements  were  small, 
their  manifestations  in  the  iterative  schemes  were  acute.  The  therrao- 
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dynamic  functions  in  the  region  of  overlapping  were  then  plotted  and  so 
adjusted  as  to  form  a  continuous  fit.  The  final  tabular  internal  energy 
values  are  based  on  the  ground  state  of  the  molecule. 

The  log10( v/Vq)  values  in  the  table  are  normalized  to  the  air 
model  of  Reference  3  at  standard  conditions.  By  changing  the  initial 
density  constant  in  sub-routine  IB,  the  normalization  of  the  table  is 
effectively  changed  to  that  of  any  initial  density  within  the  limits  of 
the  tabular  data. 

The  sub-routine  functions  as  follows  (see  Figure  l) :  Given 
a  value  of  v/vo  (normalized  to  the  table)  and  E,  the  sub-routine  searches 
the  table  values  of  10Siq(v/Vq)  and  finds  the  three  nearest  ones.  This 
effectively  locates  the  three  bracketing  sets  of  PV,  E  values  of  constant 
log10(v/V0)  or,  let  us  say,  curves  1,  2,  and  3.  Having  found  these  bracket 
ing  curves,  it  then  searches  the  sets  of  PV,  E  values  of  each  bracketing 
curve  for  the  three  PV,  E  sets  nearest  the  given  value  of  E.  It  then 
quadratic ally  interpolates,  for  each  curve,  the  three  PV,  E  values  for 
the  value  of  PV  corresponding  to  the  given  E.  Finally,  it  quadratically 
interpolates  the  three  above  determined  PV  values  and  their  corresponding 
l°glC)( V/V0)  values  for  the  PV  corresponding  to  the  given  V/VQ.  The 
pressure  is  then  obtained  from  this  PV  value,  since  V  is  known. 

10.  INITIAL  CONDITIONS  (STARTING  DATA) 

The  initial  conditions  for  this  integration  were  chosen  to  be 
those  of  Problem  M  (see  LA-2000,  Chapter  6,  Section  6.2).  These  condi¬ 
tions  were  prepared  by  J.  Hirschf elder  and  J.  Magee  and  represent  the 
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Given  E 


Schematic  diagram  of  double  quadratic  interpolation 
X  nearest  points  on  bracketing  curves 

0  interpolated  PV  values  corresponding  to  given  E 
for  each  bracketing  curve 

&  interpolated  value  corresponding  to  given  E  and 
given  V/Vn 


hydrodynamic  state  of  a  15*5  KT  explosion  12  milliseconds  after  detonation. 
Radiation  transport  was  considered  in  the  determination  of  the  fire  ball 
developmental  stages.  The  pressures,  velocities,  radii,  and  volumes  of 
Problem  M  were  used  directly;  however,  the  internal  energies  were  de¬ 
termined  from  the  present  equation  of  state  using  the  given  pressures  and . 
volumes.  This  was  necessary  in  order  to  maintain  a  pressure,  a  volume, 
and  an  internal  energy  consistent  with  this  equation  of  state.  Similarly 
the  shock  pressure  of  Problem  M  was  used  directly  and,  from  the  Hugoniot 
conditions,  using  this  shock  pressure,  the  shock  volume  and  energy  were 
determined . 

These  initial  data  are: 

Undisturbed  Medium 

p.  =  1.1613  x  10-5  gtn/cj,,?  p  =  1  x  10^  dynes/cm2  En  =  2.1421  x  10^ 

ergs/gm 

Shock  Quantities 

ry  O  1  0  I 

\|r  =  7.7246  x  10'  dynes/cm  =  I.3860  E  =  3-1162  x  10  ergs/gm 

4  =  7.9879  x  10  meters  4  =  2 .7608  meter s/millisec 

V|  =  2.7253  x  10" 5/sec 
Regular  Mass  Points 

Initially  there  were  16  mass  points.  The  isothermal  sphere 
is  contained  between  the  center  and  mass  point  1. 

-2 

Ar  =  3*994  meters  At  =  6.25  x  10  millisec 
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7.4391  X  10'  3-5017  X  10  7.9879  X  10  2.3823  X  10  1.5558  X  10' 


These  data  are  identical  to  Problem  M  starting  data  except  for 
the  change  to  a  new  equation  of  state,  which  changes  the  total  energy  to 
13.2  KT;  also.  At  is  one-fourth  that  of  Problem  M,  and  there  are  cor¬ 
responding  changes  in  the  velocity  at  n-^. 

11.  RESULTS  AND  COMPARISON  TO  PROBLEM  M 

The  re suits  of  this  integration  are  presented  in  graphical  form 
in  Figures  2-16.  Comparison  to  Problem  M  is  also  shown.  The  integration  of 
the  present  problem  was  carried  out  to  a  shock  overpressure  of  about  3 
atmospheres,  corresponding  approximately  to  a  shock  radius  of  350  meters 
and  an  elapsed  calculation  time  of  300  milliseconds  or  312  milliseconds 
after  detonation.  The  calculation  was  not  carried  further,  however,  as 
the  large  amount  of  machine  time  required  indicated  a  need  for  a  logis¬ 
tical  change  in  this  problem.  The  time  difficulty  arose  not  through  the 
use  of  the  tabular  equation  of  state,  but  through  the  ever  increasing 
number  of  mass  points  added  to  the  calculation  as  the  shock  progressed. 

It  is  believed  that  if  the  methods  of  H.  H.  Goldstine  and  J.  von  Neumann^ 
were  employed,  this  difficulty  would  be  alleviated. 

The  comparison  to  Problem  M  is  quite  good,  as  the  comparative 
figures  indicate.  The  small  differences  in  the  comparative  plots  are 
due  primarily  to  the  differences  in  the  equations  of  state  used.  Also, 
the  method  of  adding  new  fluid  elements  contributes  to  some  of  this 
difference  (see  Figure  1 6).  The  number  heading  each  adiabat  in  Figure  1 6 
refers  to  the  fluid  element  followed.  Initially,  since  there  are  1 6  fluid 
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elements,  the  pres sure -volume  point  heading  the  adiabats  are  identical; 
however,  the  pre s sure -volume  points  of  the  newly  added  points  (17,  etc.) 
do  not  agree  because  of  the  difference  in  the  adding  schemes. 

Since  it  would  be  extremely  tedious  to  compute  the  energy  of 
Problem  M  for  each  integration  cycle,  no  comparative  plot  is  made  in 
Figures  13,  14,  and  15.  The  initial  nonsmoothness  of  the  internal  energy 
curve  can  be  attributed  to  the  two  point  interpolation  scheme  used  for  the 
energy  calculation  in  the  region  between  the  shock  and  the  last  mass  point. 
Although  the  starting  data  of  Problem  M  had  an  energy  content  of  13.5  KT, 
that  of  the  present  problem  contained  13 . 2  KT.  Again  this  difference  is 
due  to  the  equation  of  state  differences.  The  positive  change  in  the 
total  energy  from  the  first  integration  cycle  to  the  last,  representing 
about  1400  integration  cycles,  was  0.66  per  cent.  The  total  energy  of 
Problem  M,  with  the  shock  front  located  at  2000  meters,  is  13*1  KT, 
representing  a  3  per  cent  loss. 
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Fig.  2.  Pressure  as  a  function  of  the  Eulerian  radius  at  t  =  0  milliseconds  for  both  Problem 
and  Problem  Herman 


Pressure  as  a  function  of  the  Eulerian  radius  at  t  =  10  milliseconds 
•  Problem  M,  °  Problem  Herman 


R  (meters) 


Fig.  4.  Pressure  as  a  function  of  the  Eulerian  radius  at  t  :  30  milli¬ 
seconds;  •  Problem  M,  0  Problem  Herman 
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Pressure  as  a  function  of  the  Euler ian  radius  at  t  =  50  milliseconds;  •  Problem 
o  Problem  Herman 


Pressure  as  a  function  of  the  Eulerian  radius  at  t  =  100  milliseconds;  •  Problem 
°  Problem  Herman 
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Fig.  11.  Shock  pressure  as  a  function  of  time;  -  Problem  M, - Problem  Herman 
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Fig.  12.  Shock  velocity  as  a  function  of  time;  - Problem  M, - Problem  Herman 


(8I0I  X  S&ia) 


k6 


Fig.  15*  Total  kinetic  energy  as  a  function  of  time 


Hugoniot  curve  and  some  adiabatic  curves  in  the  pressure-volume  plane 
- Problem  M, - Problem  Herman 


